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Abstract 

Based on a recently developed variational method, we explore the proper- 
ties of the Holstein polaron on an infinite lattice in D dimensions, where 
1 < D < 4. The computational method converges as a power law, so that 
highly accurate results can be achieved with modest resources. We present 
the most accurate ground state energy (with no small parameter) to date 
for polaron problems, 21 digits for the one-dimensional (ID) polaron at in- 
termediate coupling. The dimensionality effects on polaron band dispersion, 
effective mass, and electron-phonon (el-ph) correlation functions are investi- 
gated in all coupling regimes. It is found that the crossover to large effective 
mass of the higher-dimensional polaron is much sharper than the ID polaron. 
The correlation length between the electron and phonons decreases signifi- 
cantly as the dimension increases. Our results compare favorably with those 
of quantum Monte Carlo, dynamical mean-field theory, density-matrix renor- 
malization group, and the Toyozawa variational method. We demonstrate 
that the Toyozawa wavefunction is qualitatively correct for the ground state 
energy and the 2-point electron-phonon correlation functions, but fails for the 
3-point functions. Based on this finding, we propose an improved Toyozawa 
variational wavefunction. 

PACS numbers: 74.20.Mn, 71.38.-hi, 74.25.Kc 
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I. INTRODUCTION 



The Holstein model, as a paradigm for polaron formation, has attracted renewed inter- 
est in recent years because several lines of experimental evidence support the presence of 
polaron carriers in strongly correlated electronic materials, including colossal magnetoresis- 
tance manganites [0, organics 0, quasi-ld systems, and high-Tc cuprates [0,3. Theoretical 
research on polaron physics began six decades ago, and the problem remains unsolved due to 
its intrinsic many-body complexity from the electron-phonon interaction. (The problem of 
excitons coupled to phonons is formally equivalent 0.) Standard perturbation treatments 
I^J^ are usually limited to a particular parameter regime. With constantly growing com- 
putational resources, various numerical techniques have been applied to polaron problems 
in recent years, which give the most reliable results in the physically interesting crossover 
regime. These techniques include finite-cluster exact diagonalization (ED) P^I^], quantum 
Monte Carlo (QMC) density-matrix renormalization group (DMRG) |]l6l, and the 

global- local variational method (GLVM) ||17|| . 

Recent numerical studies have focused on the ID lattice model. Due to the enormous 
phonon Hilbert space in three dimensions, the dimensionality effects on the polaron prob- 
lems are less studied except in the adiabatic (or semiclassical) approximations |jl8|,|19[, and 



in perturbation theory |^|. QMC is also capable of computing the energy and effective 
mass of the 3D polaron, but the full dispersion E{k) is only reliable in the strong-coupling 
regime. However, with a recently developed variational method, we can compute the polaron 
effective mass, band dispersion, and el-ph correlation functions of the ground and low-lying 
excited states in all coupling regimes, preserving the full quantum dynamical feature of 
phonons (details can be found in Ref. [^). The variational space is defined on an infinite 
lattice, although only a finite separation is allowed between the electron and the surround- 
ing phonons in current implementations. We systematically expand the variational space so 
that highly accurate results can be achieved with modest computational resources. 

The main purpose of this paper is to characterize the Holstein polaron in higher dimen- 
sions. We consider a single-electron Holstein Hamiltonian on a D-dimensional hypercubic 
lattice, 

H = He + H el-ph + Hph 

= — t ^ (cjcj + h.c.) — A ^ cjcj(aj + aj) + u; ^ a]aj , (1) 
<*j> i j 

where c] creates an electron and a] creates a phonon on site j. The parameters of the model 
are the nearest-neighbor hopping integral t, the el-ph coupling strength A, and the phonon 



frequency uj. The electron is coupled locally to a dispersionless optical phonon mode [22 
There are two commonly defined dimensionless control parameters, the adiabaticity ratio 
7 = uj/t, and the el-ph coupling strength a = Ep/2Dt, which is defined as the ratio of 
polaron energy for an electron confined to a single site Ep = A^/cu, and the free electron half 
bandwidth 2Dt. The strong (weak) coupling regime refers to a > 1 (< 1), and the adiabatic 
(antiadiabatic) regime refers to 7 < 1 (> 1). An additional dimensionless parameter is 
g = X/uj, which appears in strong-coupling perturbation theory. 

A variational space is constructed beginning with a root state, the electron at the origin 
with no phonon excitations, and acting repeatedly with the off-diagonal terms {t and A) in 
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the Hamiltonian, Eq. [T|. States in generation m are those that can be created by acting 
m times with off-diagonal terms. All translations of these states on an infinite lattice are 
included, and the problem is diagonalized for a given momentum k using a Lanczos method 
pT| . A small variational space, with 7 states per electron site (unit cell), is shown in Fig. 
|1]. (The more accurate numerical computations are done with over 10'' states per unit cell.) 

The total number of states Ngt per unit cell after m generations increases exponentially, 
as {D + 1)™, where D is the spatial dimension. (The bipolaron has the same {D + I)"* 
dependence, but with a larger prefactor.) The perhaps surprising fact is that while the 
size of the Hilbert space grows exponentially with m, the error in the ground state energy 
decreases exponentially, because states are added in a fairly efficient order. Figure ^ shows 
the fractional error in the ground state energy as a function of the number of basis states 
in Hilbert space. The accuracy is determined by comparing the energy as the size of the 
Hilbert space is increased. At intermediate coupling in any dimension, the energy improves 
by about a factor of 8 with each generation In ID, each added generation approximately 
doubles the size of the Hilbert space, whereas in 4D, the size increases five-fold. This rapid 
convergence at intermediate-coupling is valuable since no analytic approach is reliable in this 
regime. Table I lists the energies for ID to 4D polarons at intermediate- to weak-coupling. 
The accuracy, 21 digits for ID polaron, is high compared to that of other numerical methods, 
such as 2 or 3 digits for QMC, 6 digits for DMRG (or GLVM), and up to 8 digits for ED 



2J]. Moreover, for the 3D polaron at intermediate- to strong-coupling, an energy accuracy 
of 8-10 decimal places can be achieved in the nonadiabatic regime with fewer than 3 x 10^ 
basis states. To obtain an accuracy beyond 13 digits, the code is executed in quadruple 
precision. The present variational method requires only power-law time to achieve a given 
accuracy (in any dimension), which is a qualitative improvement on exact diagonalization 
as it is currently implemented, the latter requiring exponential time. 

In this paper we present detailed studies of the dimensionality effect on the Holstein 
polaron. First of all, we explore the polaron characteristics in the k = ground state, and 
compare our results with previous studies from QMC, DMRG, and dynamical mean-field 
theory (DMFT). Secondly, we compute the el-ph correlation function and the polaron energy 
E{k). Finally, the validity of the Toyozawa variational method is investigated by calculating 
the ground-state energy, and the 2-point and 3-point el-ph correlation functions. 



II. SMALL-POLARON CROSSOVER 

A. Quasiparticle weight Zj. and effective mass m* 

The small polaron crossover or "self-trapping transition" has been one of the core issues 
in polaron problems. Adiabatic theory suggests that the polaron in 2D and 3D (but not 
in ID), is in an "extended" state with an infinite radius below an el-ph coupling threshold 
Ac, and beyond which is a "localized" state with infinite effective mass. (This phenomenon 
is usually termed the "self-trapping transition".) However, our studies confirm that in 
all dimensions, there is a crossover rather than a self-trapping transition (ground state 
properties are analytic), if the parameters are finite. This result is consistent with some 
other recent studies [p0| , p!4[ and corroborates the theorem of Gerlach and Lowen p5 . 
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The quasiparticle weight (renormahzation factor) is defined by the overlap (squared) 
between the bare electron and a polaron, i.e., 



Zfc = |(v[/o,fc|4|0)f , (2) 

in which |^o,fc) is the ground state wavefunction of a polaron and |0) is the vacuum state. 
can be measured in photoemission or tunneling experiments. Figure ^ shows the crossover 
of Zj:^Q as a function of 1/a at g = 3, for ID to 4D cases. We see that the crossover to 
large effective mass of the higher- dimensional polaron is much sharper than the ID polaron. 
For D > 1, a. fairly abrupt crossover occurs at a > 1, whereas the crossover for the ID 
polaron spans a wide range of a. With a smaller g (but greater than 1), the crossover 
will be slower but with the same dimensional characteristics. In the limit 1/a —>■ 0, the 
phonon wavefunction contracts to the electron site, with = exp(— f/^). The inset shows a 
comparison of Zj. and mo/m* for the ID polaron. Their fractional difference S, defined as 
(mo/m* — Zk)/Zk, is shown as a dotted line. The maximum 6 is 22%, in the intermediate 
coupling regime, while the minimum occurs as 1/a ^ (small t), where S is the order of 
from strong-coupling perturbation theory (SCPT). We find that 6 decreases significantly as 
the dimension increases. The maximum difference Smax is 4.5% for the 2D, and 2.0% for the 
3D polaron. For g = a/5, 6max in 3D drops to 0.63%. 

The ground state energy E satisfies E = — 2tcos(fc) + J]{k,E), where J]{k,E) is the 
self-energy. Z^^q is the probability of the wavefunction on the root site, and from first order 
perturbation theory Zk=o = dE/deo, resulting in Zk=o = l/[l — dT,{k, E)/dE]. The origin of 
the difference between the inverse mass and Zk=o lies in the ^-dependence of the self-energy, 

mo ^ _ 1 d'nkE) enk^ 

where the derivatives are evaluated at the ground state energy E = Eq and k = 0. In the 
variational space of Fig. |I] or in the full space, the self-energy has nonzero ^-dependence 
because distinct unit cells are connected at branch level (path 1 — 2 — 3 . . .), in addition to 
the trivial connection at root level. A restricted variational space, the comb basis, allows 
phonon excitations only on the electron site, as shown in Fig. |[ In this subspace, the self- 
energy is k-independent, since the only path between unit cells is at the root level. The 
self-energy remains k-independent even in a larger space in which the tree trunks sprout 
lateral branches, so long as the branches do not connect to neighboring unit cells. For these 
cases, the Z factor and inverse mass are identical, 6 = 0. In 0{t) SCPT, S vanishes for the 
same reason and Zk=o = mo/m* = exp(— (7^). 

The effect of dimensionality on 6 is made plausible by the following. In the Holstein 
model, the dimensionality D does not directly affect the term Hgi.ph in Eq. 1, because the 
el-ph coupling is local and the phonon is dispersionless. High dimensional polarons share the 
same simplicity of the el-ph coupling as ID. Furthermore, we see (in next section) that the 
el-ph correlation length decreases as D increases. Thus the k-dependence of the self-energy 
weakens in higher dimensions. The above arguments do not, however, hold for the Frohlich 
model (or the extended Holstein model) with longer-range el-ph coupling pB|-P^, where Zk 
and mo/m* behave quite differently. 
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B. Comparison with QMC, DMRG and DMFT 



Figure |^ shows our results for the effective mass as a function of a (at fixed 7 = 1.0) in 
comparison with DMRG and QMC. Our results are accurate to at least four digits, which is 
well below the linewidth. In all cases, Fig. |^(a)-(c), m* /uiq increases slowly when a is small, 
followed by a rapid increase when a > 1. Since it is calculated at = t = 1.0 (not a small 
t), the mass behaves differently than exp((7^) that SCPT suggests. Note that the crossover 
is more rapid as D increases, which is consistent with the results in previous section. In 
every dimension, our results are in quantitative agreement with QMC. The numerical error 
in QMC is approximately 0.1% to 0.3% |M, which is good though less accurate than finite 
cluster ED or the present approach. DMRG is generally considered a powerful tool in 
dealing with many-body problems. Using DMRG, Jeckelmann and White have calculated 
Holstein polaron properties in ID and 2D. DMRG seems to be most successful calculating 
the ground state energy (at = 0) and el-ph correlation functions. However, finite-size 
scaling is required for DMRG to compute m* ||T6[, which becomes more difficult for D > 1. 
In ID, Fig. 0(a), the results from DMRG are as accurate as QMC. DMRG does not, however, 
calculate the mass accurately in 2D, Fig. |^(b). 

Dynamical mean-field theory has previously been applied to the Holstein polaron problem 
||29|| . The approach is exact in infinite dimensions but an interpolation to 3D lattices is 
made possible by using a semielliptical free density of states N{E) to mimic the low-energy 
features. Figure |^ shows a comparison of our results on a cubic lattice to DMFT, which is 
made by setting the bandwidths equal. Overall, in panel (a), we see a qualitative agreement 
between the two calculations. DMFT is accurate in the strong-coupling regime, where the 
surrounding phonons are predominately on the electron site. This is also the regime where 
strong-coupling perturbation theory works well. In Fig. ^(b), we see our numerical results in 
agreement with weak-coupling perturbation theory in A. However, DMFT fails to compute 
m* correctly in the weak-coupling regime. The reason is that in DMFT, the lattice problem 



is mapped onto a self-consistent local impurity model p0|j31[1, which preserves the interplay 



of the electron and the phonons only at the local level. We will see that the spatial extent of 
the el-ph correlations increases as the el-ph coupling decreases, which explains the significant 
discrepancy in the weak-coupling regime. It is also worth noting that DMFT neglects the k 
dependence of self-energy, i.e., the inverse effective mass is always equal to the quasiparticle 
weight. As we have pointed out above, the difference between mo/m* and Zk is not negligible 
in the intermediate- to weak-coupling regime. 



III. ELECTRON-PHONON CORRELATIONS 

Next, we compute the correlation function between the electron and the phonon dis- 
placement (lattice deformation) in the ground state, 

X{i-j) = {"^olchiaj + a])\^o). (4) 



This correlation function can be considered as a measure of the polaron size |^2[. It should 
not be confused with the "polaron radius" in the extreme adiabatic limit, which refers to 
the spatial extent of a symmetry-breaking localized state. We would like to emphasize that 
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a comprehensive study of el-ph correlation for the ground state of the 3D polaron has not 
yet been reported by any other modern numerical technique |^3|. The on-site correlation 
has been studied by DMFT and the results are compared in Fig. |^ [Q]. The on-site lattice 
distortion x(0) is shown as function of a and the rest of parameters are the same as Fig. |^. 
In Fig. 1^, x(0) is normalized to 1 when a is infinite (i.e. t ^ 0) according to limx(O) = 2g. 

Again, we obtain good qualitative agreement. The curves show an abrupt change in slope 
only for large g, where the discrepancy with DMFT is largest. 

Figure ^ shows the effect of dimensionality on the correlation function x(i — j). In the 
strong-coupling regime, panel (a) shows, in every dimension, a sharp drop on the first two 
sites and an exponentially decaying tail. For the 3D polaron at a distance of 3 lattice sites, 
x(3)/x(0) drops below 10^^. In the weak-coupling regime, panel (b), x has nearly a simple 
exponential decay with a less steep slope, which implies a nontrivial extent of the el-ph 
interplay in space. In both panels, we observe a common trend that x decays more rapidly 
as the lattice dimension increases, i.e., the surrounding phonons are more localized near the 
electron in higher dimensions. This feature enables DMFT to give sensible results in higher 
finite dimensions. We have also investigated other 2-point el-ph correlation functions such 
as {clciOjaj) (not shown), which has dimensional characteristics similar to x- 

The rapid decay of the el-ph correlation function for the higher-dimensional polaron 
suggests that the off-site el-ph interplay is relatively weak in large D. One would then 
expect the comb basis of Fig. ^, a subspace of the full Hilbert space, to give a better 
approximation in large D. We check this assumption by numerically calculating the fraction 
of the probability density in the exact ground state that resides in the comb subspace, 

Pcomb= i^olPl'^o), (5) 

where P is the projection operator onto the comb subspace and the wavefunction \E'o is 
obtained in the full variational space. Figure ^ shows Pcomb as a function of the inverse bare 
coupling constant 1/a for 1-4D cases. In both of the limits a = or a = oo, Pcomb goes to 
1. The minimum overlap occurs in the crossover regime. As expected, Pcomb gets closer to 1 
as D increases. For the 3D polaron, the minimum of Pcomb is 91.1%, in contrast to 45.8% for 
ID. These trends can also be seen analytically. In the adiabatic limit (tu = 0), perturbing in 
t from a self-trapped state with energy Ep, the self-trapping transition occurs at a = 1 — 
The leading order correction of Pcomb for the self-trapped polaron state is 

^comb = 1 Pcomb 

In the non-adiabatic limit, Acomb can be calculated by SCPT to second order in the hopping 
t. It takes the following form: 

^4g-2g2 oo oo g2{n+m) ^ 
^comb — TTn, 2~ i i~ 7 , ^2' v"^) 

nl ml [n + my 

The above expressions show that for a given a and the discrepancy Acomb decreases as 
D increases and eventually vanishes in infinite D. The comb basis should thus give a good 
account for the Holstein problem in large D. We see in Fig. |^, however, that dimension 3 is 
not high enough for the comb to give quantitatively accurate results, and that dimension 4 
is not much better. Convergence to higher dimensions is slow. 
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IV. ENERGY DISPERSION E{K) 



Most of the recently developed numerical methods are capable of computing the polaron 
band dispersion in ID. For the 2D polaron, the only non-perturbative calculations of band 
dispersion published so far were computed by finite-cluster ED |11| and Quantum Monte 



Carlo [T^. Due to the huge phonon Hilbert space in high dimensions, the previous ED 
results are limited to small clusters, so that the band dispersion can only be evaluated at a 
few k points. The QMC allows calculation of energy at any desired k point, but is limited 
to the condition that the polaron bandwidth is much smaller than the phonon frequency, 
which corresponds to the strong-coupling regime. 

The present variational approach, however, is not subject to any of the above restrictions 
pT| . Figure |iy(a) shows the evolution of the band dispersion for the 3D polaron along 
symmetry directions in the Brillouin zone at various el-ph coupling constants A. Figure 
p!0|(b) shows the corresponding Zk- Starting with weak coupling A = 2.0 (dashed line), 
the polaron band is close to the bare electron band at lower band edge. The deviation 
between them increases as k increases. When E{k) — E{0) approaches u, we observe a 
band flattening effect, similar to the ID and 2D cases, accompanied by a sharp drop of 
quasiparticle weight Z^. The large k lowest energy state can be considered roughly as "a 
k = polaron ground state" plus "an itinerant (or weakly-bound) phonon with momentum 
fc". It is the phonon that carries the momentum so as to make Zk essentially vanish and 
give a bandwidth -E(vr) — -E'(O) = u. Due to the large extent of the el-ph correlations in 
the flattened band, our results are less accurate in the flattened regime ||3^. In the case 



of intermediate coupling A = 3.5, the polaron bandwidth is narrower than the phonon 
frequency. The upper part of the band has much less dispersion than the lower part but 
with a substantial Zk . This indicates a distinct mechanism for the crossover as a function 
of k. In the case of A = 4.5, the strong el-ph interaction leads to the well-known polaron 
band collapse and a significant suppression of Zk at all k. 

V. TOYOZAWA VARIATIONAL METHOD 

Four decades ago, a simple and intuitive variational approach to the ID polaron problem 
was proposed by Toyozawa [^. This method has been successfully applied to various fields 



and revisited in a number of guises [0,0 throughout the years. It is generally beheved to 
provide a qualitatively correct description of the polaron ground state, aside from predicting 
a spurious discontinuous change in the mass at intermediate coupling. We show below that 
although the Toyozawa wavefunction gives a good account of the ground state energy and 
the 2-point functions, it fails to correctly describe the 3-point functions. 
The Toyozawa wavefunction is written as a product of coherent states, 

i^T(fc))=Ee^'''4io) niw , (8) 

j m 

where \zi) is a coherent state of the phonon mode on site i. In the antiadiabatic limit uj/t 
oo, this wavefunction gives the exact solution c]|0)|zj), where Zj = X/ou and the other z's are 
zero. For the general case, momentum k = 0, the z's are real and symmetric: Zj^m = ^j-m- 
To determine the validity of the Toyozawa wavefunction, we probe the structure of the 
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phonon cloud in the k = ground state by computing the following 2-point and 3-point 
el-ph correlation functions, 



a2(j) = {clcoa]aj) , 
Oi^{i,m) = {clcoa!jajal^arn)- 



(9) 
(10) 



The z's in Eq. |^ are optimized so as to give a minimum energy. It can be proved that the 
optimal z's decay exponentially as a function of el-ph separation. Thus it always gives purely 
exponentially decaying 2-point functions regardless of the el-ph coupling. This, however, is 
not true of the numerically exact results. Table II and Fig. compare Toyozawa's and 
the numerically exact results for intermediate coupling B^. We notice that the Toyozawa 



wavefunction gives reasonably accurate results for the ground state energy, 2-point functions, 
and Zk=o. In Table II, the fractional error in energy is about 1% (with Zj-^i/zj = 0.35568 
and Zfj = 0.57033). However, it gives wildly inaccurate 3-points functions. For example, the 
Toyozawa a3(l, —2) is a factor of 36 too large and 0:3(1,2) is a factor of 2 too small. The 
Toyozawa 03(5, —6) is too large by 6 orders of magnitude. This failure indicates that the 
electron does not organize its surrounding phonon cloud in the way that Toyozawa suggested. 
Instead, by directly analyzing the exact ground state wavefunction, we find that the electron 
organizes its surrounding phonons like a traveling salesman does, namely, the polaron favors 
the phonon configuration with a shorter creation path. (The length of the creation path 
is the number of off-diagonal operations, phonon creations and electron hops, required to 
create a state from the bare electron state. Shorter paths are favored at intermediate or 
weak el-ph coupling, although more on-site phonons can be favored at large coupling.) For 
example, we have. 
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in the 2-phonon subspace. The amplitude attenuates rapidly as the phonon-creation path 
increases. 

The numerically exact result in Fig. 0(b) shows that it is far more favorable to create 
two phonon excitations on the same side of the electron than on opposite sides. Therefore, 
we propose to write a polaron as a sum of two asymmetric clouds, one extending like a 
comet-tail primarily off to the right and the other extending to the left. 



\%{k)) = Ee^'^cJlO) (. . . 1^,-2) k,-i)k,)k,+i)k,+2) . . . + 
j 

■ ■ ■\Zj+2)\Zj+l)\Zj)\Zj-l)\Zj-2) ■ ■ ■) , 



(11) 



where Zj-m 7^ zj+m, and the normalization factor has been omitted. The optimized (min- 
imum energy) phonon wavefunction in Eq. |ll| is strongly asymmetric, and in fact changes 
sign on one side, as shown in Table III. The main purpose of Eq. O is to investigate how 



the simplest asymmetric wavefunction improves the Toyozawa method. Shore and Sander 
have proposed a more complicated wavefunction 1^1^55) which is a sum of the symmetric 
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term in Eq. |] and the two asymmetric terms in Eq. |TT] . (Asymmetric wavefunctions 
are also considered in Ref. |^|.) The number of independent variables in and 



is ^A^, A^, and |A^ respectively, where A^ is the number of sites that allow phonon exci- 
tations. The minimum energies from the above methods are compared in Fig. [l^. It is 
clear that the energies are improved as we expand the variational space, C C ^I^s'^. 
The Shore-Sander method shows the most substantial improvement in the crossover regime. 
The comparison of other polaron properties is shown in Table II and Fig. |TT[ Our trial 
wavefunction improves the energy by 30% and the k = Z-factor by 50% compared 
to 75% and 66% respectively from Shore-Sander wavefunction. In Fig. |TT](a), Eq. ^ gives a 
more accurate 2-point function 02 (j) than the original Toyozawa. It similarly improves the 
other 2-point function Xj (not shown). Panel (b) shows two 3-point functions, asiJjj + 1) 
and 03 (j, —j — 1). Due to its symmetric phonon cloud, the Toyozawa wavefunction must 
give exactly the same result for the two 3-point functions. In contrast, the exact results 
show that a3(j, j + 1) ^ «3(j, —j — !)• Equation ^ gives the correct behavior of the two 
3-point functions on nearby sites, but loses quantitative accuracy in the tails. Although 
the Shore-Sander energy is better than that of Eq. [11], the Shore-Sander 3-point functions 
are actually worse. The simplest attempt {"^'t) to correct the identified shortcomings in the 
Toyozawa variational wavefunction appears to be a step in the right direction, although not 
as quantitatively accurate for most properties as variational methods with more parameters 

mm. 



VI. CONCLUSION 

In summary, we have performed extensive numerical studies of the Holstein polaron in 
spatial dimensions 1 through 4. The numerical method used adds basis states to the Hilbert 
space in an efficient order, resulting in an error that scales as a power of the size of the 
Hilbert space N~f^ , where ^ is a nonuniversal exponent 3 at intermediate coupling in ID, 
and ~ 1.6 in 3D. This is a qualitative improvement over standard exact diagonalization, 
which requires exponential effort to achieve a given accuracy. Using modest computational 
resources, we obtain by far the most accurate polaron energies and wavefunctions available 
from ID to 4D at intermediate coupling. 

Previously, a thorough investigation of the dimensionality effect, including correlation 
functions, was out of reach of numerical methods. The main findings of the dimensionality 
effects on the the Holstein polaron are summarized as follows: The crossover from quasi-free 
to large effective mass is found to be much sharper in higher dimensions. As was recognized 
previously, there is no symmetry-breaking self-trapping transition for finite parameters in 
any dimension, as suggested by adiabatic theory (although there is a phase transition in the 



first excited state |]2l[). See also Ref. Our results for m* agree with QMC, although there 
is a discrepancy with DMRG in D > 1. The electron-phonon correlation functions decay 
significantly faster in higher than lower dimensions. This implies a shorter el-ph correlation 
length in large dimensions and leads to a diminishing difference between the inverse effective 
mass mo/m* and the wave function renormalization Zj:^^ as D increases. The DMFT 
approach thus gives better results in higher dimensions. Our comparison shows that DMFT 
gives qualitatively correct results for the effective mass, mean phonon number, and on-site 
phonon distortion in the intermediate- to strong-coupling regime. We also examine the 
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comb-basis approach which hmits the el-ph correlation to the on-site level as DMFT does. 
The discrepancy between the comb basis and the full basis decreases slowly as D increases. 

Finally, our approach is compared to the well-known Toyozawa variational method. We 
quantitatively examine the method in the intermediate-coupling regime. Overall, the Toy- 
ozawa wavefunction gives reasonably accurate energy and 2-point functions, but fails seri- 
ously for the 3-point functions. (The numerically exact 3-point functions are quite different 
for excitations on opposite sides of the electron compared to the same side, whereas the 
Toyozawa wavefunction predicts that they are identical.) We propose an improved varia- 
tional wavefunction, a sum of two asymmetric phonon clouds (Eq. |r^), which gives improved 
3-point functions, and somewhat more accurate results for the energy, Z-factor, and 2-point 
el-ph correlation functions. 

For all the polaron features calculated, the present numerical approach compares fa- 
vorably to other numerical methods in terms of accuracy, ease of implementation, and the 
ability to compute ground and excited state energies and correlation functions. It can also 
be directly applied to study the effects of dimensionality on other interesting problems, 
such as the Frohlich model or extended Holstein model with longer range electron-phonon 
interactions, and to bipolaron problems. 
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TABLES 





ID 


2D 


3D 


4D 


Eq 


-2.46968472393287071561 


-4.814735778337 


-7.1623948409 


-9.513174069 



TABLE I. Polaron ground state energies at A; = in ID - 4D for a = 0.5, g = 1.0, and t = 1.0 . 







"2(0) 



a3(l,2) 03(1,-2) 



-0 



this work -2.69356579774920. . . 0.40770 0.0004691 0.000005888 0.627322. 

Toyozawa -2.662819 0.32527 0.0002142 0.0002142 0.65738 

Eq. |ll| -2.671530 0.34240 0.0007649 0.000003244 0.64271 

Shore-Sander -2.685826 0.37780 0.0005572 0.0001132 0.63757 



TABLE IL A comparison of the ground state energy £'0, two- and three-point el-ph correlation 

. by the present metl 
|). Parameters are A = 1.2, u) 



functions, and Zk=o, evaluated by the present method, the Toyozawa, Eq. 11, and Shore-Sander 
wavefunctions (^'^^ in Ref. 



1, t = 1 n = 1. 
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site j 




-6 


-0.12384D-03 


-5 


-0.39019D-03 


-4 


-0.11875D-02 


-3 


-0.32178D-02 


-2 


-0.48444D-02 


-1 


0.35290D-01 





0.58515D+00 


1 


0.38153D+00 


2 


0.14043D+00 


3 


0.46112D-01 


4 


0.14632D-01 


5 


0.45908D-02 


6 


0.14349D-02 



TABLE III. A partial list of the optimized phonon wavefunction Zj in Eq. 
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FIGURES 




FIG. 1. A small variational Hilbert space, a subset of the generation 3 space, is shown for the 
ID polaron. Basis states in the many-body Hilbert space are represented by dots, and nonzero 
off-diagonal matrix elements by lines. State |1) in the root state, an electron at the origin with 
no phonon excitations. Vertical bonds create or destroy phonons. State |2) is an electron and one 
phonon, both at the origin. State |3) is an electron on site 1, and a phonon on site 0. The dots 
can also be thought of as Wannier orbitals in a one-body periodic tight-binding model. 
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FIG. 2. The fractional error A in the polaron ground state energy as a function of the number 
of basis states Ngt in the Hilbert space for parameters a = 0.5, g = 1.0, and t = 1.0. 
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FIG. 3. Quasiparticle weight -^^^q as a function of the inverse couphng strength 1/a for ID 
(sohd hne), 2D (dotted hne), 3D (dashed Hne), and 4D (long dashed Une). a is varied by changing 

the hopping t at fixed uj and A. The comb basis approximation (sec below) is shown as a dot-dashed 
line. The inset shows the comparison of the inverse effective mass mo/m* (dashed line) and Zk 
(solid line) for ID. The fractional difference S = {mo/m* — Zk)/Zi; is plotted as a dotted line. 
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l\ 14 



FIG. 4. The comb basis, a variational space in which phonon excitations are present only on 
the electron site. Vertical lines create phonons and horizontal lines are the electron hops. State 1 1) 
is an electron on site and no phonons. State |2) is an electron and one phonon, both on site 0. 
State |3) is an electron and two phonons, all on site 0. State |4) is the a translation of state 
The comb basis is a subset of the larger variational space. As in DMFT, it only keeps track of the 
on-site el-ph correlations. 
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FIG. 5. The effective mass m* for the (a) ID, (b) 2D, and (c) 3D polaron is compared to 
DMRG [|16) and QMC g calculations. (No DMRG data is available for 3D polaron.) In all cases, 
CO = 1.0, and t = 1.0. Note different horizontal scales. 
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FIG. 8. Correlation x of the electron density and the phonon displacement as a function of 
distance {i — j) for the 3D polaron along the (1, 0, 0) direction, the 2D polaron along the (1, 0) 
direction, and the ID polaron at (a) strong coupling, (b) weak coupling, a; = 1.0. Note the different 
vertical scales. 
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FIG. 9. The probability density in the ground state that resides in the comb subspace Pcomb 
as a function of the inverse bare couphng strength 1/a, for the 1-4D polaron. The parameter set 
is the same as in Fig. |^. 
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FIG. 10. Ground state energy E{k) of the 3D polaron in panel (a) and quasiparticle weight 
Zk in panel (b) for three different el-ph coupling eonstants, A = 4.5 (solid line), A = 3.5 (dotted 
line), and A = 2.0 (dashed line). Other parameters are lo = 2.0 and t = 1. The dot-dashed line 
in (a) is the dispersion of a bare electron. The corresponding ground state energies E{k = 0) are 
-10.608348, -8.0642850, and -6.588526818 respectively. 
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FIG. 11. (a) The 2-point function a2(j) is evaluated in ID by the present variational method 
(solid line with squares), the Toyozawa method (dashed line with circles), and the modified Toy- 
ozawa method Eq. |ll| (dashed line with diamonds), (b) The 3-point functions a3(j, i + 1) (solid 
lines) and a3(j, —j — 1) (dashed lines). The symbols are the same as in (a). Note that the plain 
Toyozawa method gives exactly the same results for the two 3-point functions, which in fact differ 
widely. Parameters are uj = 1.0, t = 1.0, and A = 1.2. 
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FIG. 12. A comparison of ground-state energy as a function of coupling constant from various 
variational approaches for a; = l,t = 1. 
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